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ABSTRACT 

Thermal convection phenomena In a vertically oriented Brldgnan-Stockbarger apparatus 
were modeled by computer simulations for different gravity conditions# ranging from 
earth conditions to extremely low gravity# approximate space conditions,, The model- 
ing results were obtained by the application of a state-of-the-art# transient# 
multi -dimensional » completely densl metrically coupled# discrete-element computation- 
al model which was specifically developed for the simulation of flow# temperature# 
and species concentration conditions In two-phase (solid-liquid) systems. The compu- 
tational model was applied to the simulation of the flew and the thermal conditions 
associated with the convection phenomena In a modified Germanium-Silicon charge 
enclosed In a stationary fused-slllca ampoule. The results clearly Indicated that 
the gravitational field strength Influences the characteristics of the coherent 
vortical flow patterns# Interface shape and position# maximum melt velocity, and 
Interfaclal normal temperature gradient. 


INTRODUCTION 

It is well established that melt convection exerts a strong Influence on the 
processes that control the growth of compound semiconductor crystals. In Bridgman 
crystal growth# convection Is driven by density differences In the melt which 
continuously Interact with the local gravitational field. The majority of the 
problems encountered during attempts to produce homogeneous# defect-free crystal s# 
are attributed to the direct result of melt convection through Its Influence on the 
heat and mass .transport during solidification. Consequently# measures are usually 
taken to minimize convection and to enhance the possibility of crystal growth In a 
diffusion dominant regime. The microgravity environment of space offers excellent 
controlled laboratory conditions for studying the characteristics of the convection 
phenomenon In crystal growth which Is simultaneously controlled by the density 
differences in the melt and the gravitational force field. 

The physicochemical characteristics of crystals grown In space under microgravity 
conditions are usually quite different than their earth based counterparts grown 
under otherwise similar conditions. These differences however cannot be explained on 
the basis of the currently available mathematical models with analytic solutions 
(e.g., Sukanek# 1982; Lehoczky and Szofran# 1982)# because their Inherent simplify- 
ing assumptions usually tend to discard certain Important details of the physical 
phenomena which control crystal growth. Order of magnitude studies (e.g.# Camel and 
Favler# 1984) can also yield useful Information# but again# generally lack the 
capability to consider the detailed consequences of the physical phenomena which are 
necessary to explain the different results of the experiments. Relatively simple# 
analytic or computational models which only consider diffusion (e.g., Corlell and 
Sekerka# 1979) are also usually Inadequate in studies concerned with the understand- 
ing of the physics of crystal growt 1 * •“‘oblp — ^-’use they completely disregard the 
presence of the convection phenomenor which can strongly Influence the solidification 
process. 

The overall objective Is to provide the required state-of-the-art computational 
modeling capability which can consider all the effects of the Important physical 
phenomena# Including convection and diffusion transport# body forces induced by 
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electromagnetic fields, Coriolis acceleration, and the acceleration due to transient 
motion, both under earth and microgravity space conditions. An Important part of the 
overall objective Is to understand the role of melt convection In crystal growth. 
Considering the fact that melt convection has not been established to be universally 
detrimental, a high level understanding of the physical phenomena can possibly 
provide a means of controlling convection, and hence, optimizing or Improving the 
crystal growth process. 


This study presents the results of the preliminary applications of the computational 
model to the simulation of crystal growth In a stationary cylindrical ampoule, 
under axially symmetric, thermally driven flow conditions, without any solutal 
convection, and with the gravitational force aligned with the vertical ampoule axis. 
The results are presented for sytematlcally varying gravitational field strength 
conditions, and for typical furnace conditions, In order to Illustrate the effects 
of the gravity field on the crystal growth process In a relatively well controlled 
Qrldgpnan system. 
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SUMMARY OF DISCRETE-ELEMENT COMPUTATIONAL MODEL 

The discrete-element method (DEM) (e.g., Eraslan, 1974; Eraslan, Lin and Sharp, 

1982; Eraslan and Lin, 1985; Eraslan 1985) directly utilizes the 1 ntegral -forms of 
the physical principles In the development of the spatially discretized computation- 
al systems for the solutions of the velocity, temperature and species mass-fraction 
(concentration) conditions, by Identically satisfying the conservation principle for 
melt mass In each and every "element". Since, It does not require, any consideration 
of the continuum-limit partial differential equations (l.e., continuity equation, 
Navler-Stokos equation, energy equation, etc.), associated with the formulations of 
the conventional mathematical models, the DEM Is conceptually and mathematically 
different than the finite-difference method (FDM) and the finite-element method 
(FEM); and It Is similar (only conceptually) to classical fluld-ln-cell (FLIC) and 
Integrated-finite-difference (box) direct numerical si mul atl on methods [e.g.. 

Gentry, Martin and Daly, 1966; Harlow and Amsden, 1969; Roache, 1972). 


Composite-Space-Splitting Algorithm 

The application of the three-dimensional DEM to the simulation of the convection and 
diffusion phenomena In crystal growth considers a special discretization form. In 
cylindrical coordinates (r, 0,z), which utilizes four partially overlapping half- 
elements, In the (r,9) plane. In each element. Two half-elements are used for the 
computation of the r-component of velocity, and two half-elements are used for the 
computation of the e-component of velocity. In each element. The three-dimensional 
method calculates (1) the melt free-surface elevation In N elements of the surface 
layer, (2) the r-component of velocity In 2N half-elements, (3) the 0-component of 
velocity In 2N half-elements, (4) the z-component of velocity on N enclosure surfa-- 
ces, (5) the pressure on N enclosure surfaces, (6) the temperature In N elements, 
and (7) the species mass-fractions (concentrations) In N elements of the computa- 
tional region. The complete computational system of 7( seven) sets of spatially 
dlcretlzed equations for elements (1=1,2,. ..,1; j=l,2,...,J; *=1,2, . . .»L) can be 
represented as 


tdl scr ete- element equations for melt-surface elevation In the elements 
of the surface layer (*=L)], 

[discrete-element equations for r-component of velocity In the 2(two) 
half-elements of the Internal layers]. 
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, [discrete-element equations for 9-component of velocity In the 2(two) 
d ‘ half-elements of the Internal layers]* 

V ,, * [discrete-element equations for z-component of velocity at the 

enclosure surfaces, between the Internal layers of the elements], 


the enclosure surfaces 


* . = [discrete-element equations for pressure at t 

between the Internal layers of the elements]. 


B [discrete-element equations for temperature In the elements], 


„ [discrete-element equations for species mass-fraction (concentration) 
In the elements]. 

The discrete-element system Is formulated by considering spatially discretized forms 
of the phenomenological parts, represented as (1) "flow” for melt mass transport In 
the conservation of mass equation, (2) "conv" for "convection" transport In all 
equations, (3) "fore" for "pressure force" and "body force" (Including the "gravity 
force" and the "Coriolis force") In all momentum equations, (4) "vise" for "viscous 
stress" In all momentum equations, (5) "cond»k" for "thermal conduction" In the 
energy equation, (6) "d1ff,D" for "mass diffusion" In the species-mass-conservation 
equations, and (7) "genr" for "generation" In the energy and species-mass- 
conservation equations. 

Composite-Time-Splitting Algorithm 

Timewise numerical Integration of the DEM model Is based on a novel, explicit compu- 
tational algorithm which uses the composite-time-splitting technique. The numerical 
Integration advances the solution In time by considering (1) a single-step, two- 
time-level partial algorithm for the calculations of "conv", the convection trans- 
port terms, which guarantees the elimination of the first-order nunerlcal dispersion 
effects, and (2) a two-step, three-time-level partial algorithm for the calculations 
of "fore", the force terms, "vise", "cond,e», d1ff»D"» the nonconvectlve transport 
terms, and "genr", the generation terms. 


With the known solution set at the Initial time t , the composite-time-splitting 
algorithm of timewise nunerlcal Integration can be represented In the generalized 
difference operator notation for the phenomenological parts. 


At the Intermediate time level t 
(1) melt-surface elevation h. 


n-H=S _ ^n + 
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(2) r-component of velocity V r , 
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(4) temperature T, 
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(5) mass-fraction (concentration) C k of specles-k: 
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(7) pressure p* 
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at the solution (final) time level t n+1 = t n + At n : 
(1) melt-surface elevation h» 
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(3) 8-component of velocity V Q # 
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(6) z-component of velocity V z » 
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Phenomenological Parts of the Difference Operators 

The spatially discretized phenomenological parts of the difference operators In 
Eqs. (1) through (14) are formulated# according to the composlte-space-spl Ittlng 

algorithm of the DEM# at the two time levels# t n and t n4i *» of the composite-time- 
splitting algorithm. 

Abbreviated forms of the “flow” difference operator for the me! t- surf ace elevation 

at the initial time level t n and at the intermediate time level t , respectively, 
can be represented as: 


[a" h] n . = Eq. (15b) , with n+k replaced by n, (15a) 

L flow , J n 
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+ additional flow terms for melt mass transport]. 
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Abbreviated forms of the "conv” difference operator for the r-component of velocity 

at the Initial time level t n and at the Intermediate time level t n4 **, respectively* 
can be represented as: 


[A* V l n . , * Eq. (16b), with n+*5 replaced by n. 


(16a) 
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Abbreviated forms of the "fore" difference operator for the r-component of velocity 

at the Initial time level t n and at the intermediate time level t n+Js , respectively, 
can be represented as: 




= Eq. (17b), with n+h replaced by n. 
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Abbreviated forms of the "v1sc !fl difference operator for the r-component of velocity 

at the initial time level t n and at the Intermediate time level t n+ij , respectively, 
can be represented as: 


£ A viscVi-Jfi,j,* = Eq * ( 18b *’ With n+is re P laced n ’ 
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v ,rr’ v ,8r’ v ,zr 

The spatially discretized phenomenological parts of the difference operators In 
Eqs. (1) through (14) for (1) the 3-component of velocity# (2) the pressure (from 
the z-component of the momentum equation)# (3) the temperature# and (4) the mass- 
fraction (concentration) of specles-k are also formulated# according to the 

composite-space-splitting algorithm of the DEM# at the two time levels# t n and t # 
of the composite-time-splitting algorithm. 

The enclosure surface areas A„# A a and A and the volumes v of the elements and 

r 0 z 

half-elements# used In the formulations of the difference operators# Eqs. (16a) 
through (18b)» are formulated according to the spatial discretization of the DEM In 
cylindrical coordinates (r»0»z)» i.e.# 


A = rA8Az, A. = ArAz, S = rA6A r, v * rASArAz . ’ (19) 

r p z 

The DEM formulation of the computational model considers the mathematical system as 
completely denslmetrically coupled# in the sense that the density variations# with 
temperature and species mass-fractions# are retained exactly in all the terms of the 
governing discrete-element equations. Therefore# the DEM model does not resort to 
the use of the simplifying Bousslnesq approximation (e.g.# Eraslan, Lin and Sharp# 
1982); but It considers the variations of the density in terms of the generalized 
compressibility fraction B defined as 

P = P 0 ( 1 + f^“) = P 0 0 + 3) (20) 

V O 

such that the directional flow rates across the enclosure surfaces of the elements 
and half-elements explicitly include the local effects of compressibility# I.e., 


= rAQAZ V (1 + 8); 6 0 = ArAz V. (1 + p); G = rA8Ar V (1 + 5) (21) 

Pjt t p j t/ 0 p j Z Z 

The general abbreviated forms of the "conv" difference operator# associated with 
convective transport for the solution temperature T and species mass-fraction 
(concentration) C^# at two time levels# can be expressed as 

[A n F]? . = Eq. (22b), with n+% replaced by n, (22a) 

r A n+Js p-in _ ] rg n+i S (F n - F n ) 

^ conv rjn+J$ /, . R n+^ \ 8>r;i-*s, j,t i”N»jit ijj#^ 

i,j 5 r 

+G n+J * , (F n - F n ) 

B.rji+ij, j,A ' i+4,j,£ 


+ additional convective transport terms (F)] (22b) 
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It Is Important to point out that the special form of the convective transport term* 
formulated as the difference between the value of the flow variable at the enclosure 
surface and Its computational value In the element or half-element# represents a new 
Intensive-form (IF)# which introduces more desirable stability characteristics to 
the numerical solutions. 


Dlrectional-Transportive-Upwlnd-Interpol atlon (DTUI) 


The general computational algorithm of the discrete-element method (DEM) Incorpor- 
ates the dl rectlonal-transportlve-upwlnd-lnterpolatlon (DTUI) technique which 
represents a second-order formulation# based on (1) half-point velocity components 
at the enclosure surfaces of the elements and half-elements and (2) the specified 
time step of numerical Integration. The second-order upwind-differencing technique 
evaluates the transported values of the intensive flow properties at specific 
locations upstream of the enclosure surface as they move with the flow during the 
specified time j)er1od» l.e. # 


W S5 f 4* A V* 


***"%» J#A a t At 


(23a) 




(23 b) 


z i->s,j,a * + Az i-4.j,A = Z t- J s,j.i " V* 


(23c) 


where# 

a fc * 0.5 - 0,6 (depending on the stability characteristics). 


(245 


With the evaluated# directional -transparti ve-upwind locations# Eqs. (23a), (23b) and 
(23c)# the computational algorithm uses the multl-dlmenslonal-lnterpolatlon (MDI) 
algorithm to evaluate the dl recti onal -transport! ve-upwind values of all the fluid 
properties at the transport enclosure surface# by using the values of the fluid 
properties In all the adjacent elements and/or half-elements. 

As examples# the MDI forms of the DTUI values of the r-component of velocity and the 
temperature or species mass-fraction (concentration) can be represented as 


V r;i-4,j ,A MDI Q i-h>3 

V 


r ; t~h, j, j±l;A,A±l 


)] (25a) 


MDI [F( 
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Discussions of Convergence (Consistency# Stability# Accuracy) 

The composite-space-time-splitting# with dlrectlonal-transportlve-upwind- 
Interpol atlon and multi -dimensional -Interpol atlon# (CSTS-DTUI-MDI) algorithm tends 
to eliminate all the nonphysical contributions of the numerical dispersion (computa- 
tionally generated) effects to the transport. The application of the heuristic 
stability technique (e.g.# Roache# 1972) to the linear model systems# representing 
the Identical wave# convection and diffusion phenomena establishes that the 
convergence of the (CSTS-DTUI-MDI) algorithm to the partial differential equations 
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(.POE) Is of accuracy order At (At square) In the time step for each of the 
phenomenological terms# provided that all the required stability criteria# 

(1) Courant-Frledrlchs-Lewy (CFL)# based on the surface wave speed for maximum melt 
depth D m > (2) Courant# based on the maximum velocity Y max In the melt# and 

(3) Neumann# based on the maxlmim dlffuslvlty e fflax In the solid and/or In the melt# 


are met# l.e.# 


At D J^ 
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At- V 

max .... > i 
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It Is Important to point out that the stability of the CSTS-TDU1-MDI algorithm Is 
completely Independent of all the considerations of the characteristic nunbers# 
Including Prandtl# Schmidt# and Grashof numbers. The expanded limit of the 
discrete-element system converges directly to the standard transient three- 
dimensional PDE system, with the additional error terms that can be represented In 
the form of higher-order-derivative (HOD) terms# l.e.# 
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Each of the remaining DEM equations also converges to the appropriate PDE with At 

order accuracy In the time step. 

It Is Important to point out that In all DEM equations# the selection of the value 
of the parameter for determining the Intermediate time level t n+1 ^ 2 ' as 

Y c * 0.5 - 0.6 (depending on the stability characteristics). (29) 

results In the elimination of the Hfc-crder terms In Eqs. (27)# (28) and all other 
PDE aquations. Therefore# the CSTS-DTUD-MDI algorithm tends to eliminate all the 
numerical dispersion effects In the solution of the DEM computational model? and 
hence# It can be classified as a consistent and stable (l.e. » convergent) algorithm 
with second-order convergent accuracy In time step At. 


Consideration of the Solid and Melt Subregions and Modeling of the Phase Change at 
the Interface 

The DEM model utilizes the same computational algorithm for the calculations of the 
solutions In both the melt and solid subregions of the entire region of the charge# 
by considering two different marker values for Identifying each element as contain- 
ing either solid (marker a 0) or melt (marker a I# as liquid). The computational 
algorithm changes the marker value In an element only as a consequence of the 
phase change conditions In t’.te element# according to a kinetic phase-change 
submodel. The phase change phenomenon (solidification or melting) Is formulated by 
considering the solidification (melting) temperature and the phase diagram of the 
charge# and the local temperature and species mass-fraction (concentration) 
distributions along the solid-melt Interface# which are controlled by the overall 
heat transfer and mass transport conditions of the crystal growth experiment. The 
special formulations of the DEM computational model and the kinetic phase-change 
submodel eliminate the need to consider different numerical solutions for the solid 
and melt subregions# and more Importantly# eliminate the need to match conditions at 
the solid-melt Interface as a moving boundary. 


APPLICATION TO A TYPICAL BRIDGMAN-STOCKBARGER EXPERIMENT 

The outlined DEM computational model was applied to the simulation of a typical 
Bridgman -Stockbarger crystal growth experiment for studying the convection pheno- 
menon under varying gravitational field strength conditions. 


Description of the Apparatus 

An ampoule# Insulated on both Its top and bottom# Is placed In the middle of a 
furnace which has a large adiabatic zone (Fig. 1). The hot and cold zones are 
arranged such that parts of the ampoule wall can be maintained at approximately 
constant temperature levels above and below the charge melting point. The ampoule Is 
made of fused silica# with 0.185 cm thick walls# an outer diameter of 1.805 cm# and 
an overall length of 7.6192 cm. The stationary ampoule contains a charge similar to 
a Germanium-Silicon compound# except that the sol utal coefficient of volume 
expansion has been set equal to zero. Because the charge contracts upon solidifica- 
tion# the top of the melt Is not held In contact with the ampoule# 'but It is allowed 
to behave as a free surface# which closely simulates the actual processing 
conditions. Thermophysical properties of the system are also given In Fig. 1. 
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0.0815 cm2/s 

crystal thermal diffusivity 
5.48 g/cin3 
charge density 
1.10 (10“3) r,m2/a 
kinematic viscosity 


Fig. 1. Schematic of a Brldgman-Stockbarger System and Thermophysical Properties. 
Description of the Computer Simulations 

In preliminary applications of the general, transient, three-dimensional DEM 
Computational model to the described simulation of the Brldgman-Stockbarger crystal 
growth experiment. It was assumed that the direction of the gravitational field was 
aligned with the vertical axis of the ampoule. Furthermore, It was assumed that the 
the axial symmetry conditions could be maintained under all applied vertical gravi- 
tational force fields. 

The computer simulations were started with an arbitrary, but realistic, solid-melt 
interface location and shape In the ampoule. For a specified gravitational field 
strength, the transient computer simulations were continued until steady-state 
conditions were attained both for the location and shape of the Interface, and for 
the coherent recirculation patterns of the convection cells In the melt. 


Discussion of Results 


The gravitational field strength was allowed to vary from l.Og to 10 _8 g In four 

steps. The flow fields resulting from this systematic variation are presented In 
Figs. 2-5 In which only half of the r-z plane Is shown. The strongest circulation at 
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any gravitational level Is always In the largest cell which Is adjacent to the wall 
In the vicinity of the heater zone. Cell circulation Is counterclockwise with the 
flow upward along the wall, as denoted by the positive value of streamfunctlon <l>. 

The center of any cell Is the maximum value of ij> for that coli. Equal Increments 
are then shown between the center of each cell and the zero value on the ampoule 
wall. Adjaoent cells always rotate In opposite directions. 

As gravity decreases, both the quantatlve and qualitative nature of the flow and 
transport change. Four cells decrease to one and then Increase again to two, but the 
nature of the flow has changed In the process. The horizontally stratified cells of 
Fig. 2 give way to the unicellular flow of Fig. 4 which In turn yields to the 
vertical cell pattern of Fig. 5. This could have a substantial Influence on the 
transport process If the transport equation Is convection dominated. Semiconductor 
melts typically have Prandtl and Schmidt nunbers on the order of 0.01 and 100.0, 
respectively. Consequently, while rather vigorous convection will be necessary to 
affect heat transport, species transport Is altered at much lower melt velocities. 
Thus, while the species distribution Is nearly uniform within the cells of Fig. 2, 
species transport between cells Is by diffusion making the horizontally stratified 
field considerably different than the vertical field. 

Figure 6, a plot of the Interface for various gravity levels, shows that convection 

has no Influence on Interface shape or location at Intensities less than 10 ^g @ . At 

higher levels the’ Interface progressively moves downward Into the cold zone and 
changes shape. A recent publication (Carl son, Frlpp, and Crouch, 1984) has shown that 
changing the sign of the Interface curvature radically alters the flow field. While 
not shown. It has been calculated that the convective component of the energy 

-2 

transport Is only \S% of the total transport at l.Og^, 5% at 10 3 q* and 
approximately zero elsewhere. 

Constitutional supercooling Is important in multicomponent melts. Increasing fluid 
motion Increases the Interfaclal normal temperature gradient, as Indicated In 
Fig. 7. The bulk fluid Is more thoroughly mixed and hence a higher temperature 
gradient appears near the Interface. Again no influence Is felt at gravitational 

levels below 10”^g where conduction dominates. Fluid motion thus diminishes the 
e 

chances for constitutional supercooling to occur. 


CONCLUSIONS AND RECOMMENDATIONS 

The convection phenomenon In Brldgnann-Stockbarger experiments, under varying 
gravitational force conditions, was simulated by using a computational model,. The 
results of the preliminary applications of the computational model Indicated that, 
as expected. If gravitational force conditions can be reduced, the melt convection 
phenomenon can also be reduced significantly In the crystal growth experiment. 

However, the preliminary applications did not consider various Important effects 
which exist under reduced gravity conditions. Including (I) the effect of reduced 
melt convection on the species mass transport In the charge, (2) the effect of 
enhanced Coriolis acceleration (force), (3) the effect of enhanced surface tension, 
and (4) the effects of sudden acceleration or deceleration, that could result from 
orbital control requirements. Considering the cost of space projects. It is 
Important to take Into account all these Important effects In the design and 
planning of the crystal growth experiments for future missions. Therefore, a state- 
of-the-art, transient, three-dimensional computational model can be recommended as a 
cost-effective design tool for simulating crystal growth experiments under realistic 
and detailed microgravity conditions that would exist In actual orbital stations. 
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Fig. 6. Melt-Solid Interface Shape 
for Various Gravitational 
Field Strength Conditions. 



Fig. 7. Interf acl al Normal Temperature 

Gradient for Various Gravitational 
Field Strength Conditions. 
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